rm(list=objects())
graphics.off()
library(INLA)
library(inlabru)
library(sp)
library(spatstat)
library(ggplot2)
library(raster) # raster altitude
library(numbers) # fonction mod
library(viridis) # couleurs graphique altitude
library(gstat) # krigeage
library(dbmss) # fonction Mhat
Les données fournies sont celles de la parcelle 16 de Paracou (année 2015) auxquelles j’ai rajouté 3 variables.
J’ai également supprimé les arbres qui sortaient de la parcelle. Pour plus de précisions concernant la construction de ces données, voir le fichier “preparation_donnees_paracou.R”.
load("data/p16_avec coord_et_altitude_2015.RData")
Cette fonction crée à partir des coordonnées x et y du semis de points, tous les objets dont nous aurons besoin pour notre étude. Elle retourne le nombre d’observations (nb_obs), un SpatialPointsDataFrame de coordonnées (coord_pts), un objet ppp des coordonnées (coord_pts_ppp), un DataFrame des coordonnées (df) et un mesh (mesh) adapté aux coordonnées.
L’argument data doit être composé de deux colonnes nommées x et y. L’argument mesh = TRUE par défaut. L’argument size spécifie la taille d’un côté du domaine (par défaut = 5), et donc les valeurs max. des coordonnées entrées en argument (deux valeurs utilisées : 5 ou 500).
Sûrement du à un problème d’implémentation, nous ne sommes pas en mesure de construire un maillage régulier sur des données allant de 0 à 500 (valeurs trop élevées qui engendre un temps de calcul démesurément long). L’argument size fixé à une valeur de 500 servira donc uniquement à tracer certains graphiques en taille réelle ; si size = 500, la fonction ne produira jamais de maillage. La solution provisoire trouvée pour nos études est de considérer un domaine de taille 5x5, et de diviser les coordonnées des arbres de la parcelle par 100.
preparation = function(data, mesh = TRUE, size = 5){
sortie = list()
# spatialPointDataFrame
coord_pts = data.frame(data)
sortie$nb_obs = nrow(coord_pts)
coordinates(coord_pts) = c("x", "y")
sortie$coord_pts = coord_pts
# dataframe pour les graphiques
sortie$df = data.frame(data)
# création de d'objet ppp point pattern à partir des coordonnées
sortie$coord_pts_ppp = ppp(coord_pts$x, coord_pts$y, c(0,size), c(0,size))
if(size==5)
{ if(mesh==TRUE)
{ # mesh apadté aux données
bnd = spoly(data.frame(lon = c(0, 5, 5, 0, 0), lat = c(0, 0, 5, 5, 0)))
boundary <- list(as.inla.mesh.segment(bnd), NULL)
mesh.loc <- coord_pts
mesh <- inla.mesh.2d(loc=mesh.loc,
boundary=boundary,
max.edge=c(0.2, 0.74),
min.angle=c(30, 21),
max.n=c(48000, 16000),
max.n.strict=c(128000, 128000),
cutoff=0.008,
offset=c(0.16, 0.43))
sortie$mesh = mesh
}
}
return(sortie)}
Comme expliqué précedemment, les maillages réguliers sur des domaines dont les côtés dépassent 50 sont extrêmement longs à construire, même avec des exigences de qualité très faibles… La solution provisoire que nous allons donc adopter pour ces études est d’effectuer notre analyse sur un domaine de taille 5X5 et de diviser les coordonnées de chaque arbre par 100.
# frontière du domaine 5x5 (sens inverse des aiguilles d'une montre)
bnd = spoly(data.frame(lon = c(0, 5, 5, 0, 0), lat = c(0, 0, 5, 5, 0)))
# même frontière en objet owin
bnd_owin = owin(xrange=c(0,5), yrange=c(0,5))
# mesh régulier sur le domaine 5x5
mesh_r = inla.mesh.2d(boundary=bnd, max.edge = 0.2)
ggplot() + gg(mesh_r) + coord_fixed() + ggtitle("Mesh régulier")
# grille utilisée pour les prédictions sur le domaine
pix = pixels(mesh_r, nx = 200, ny = 200)
# coordonnées des noeuds de cette grille dans un dataframe
pix_df = as.data.frame(pix)
# fonction qui définit le min et le max de l'échelle de couleurs pour l'affichage
colsc <- function(max, min=0) {scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"RdYlBu")), limits = range(c(min, max), na.rm=TRUE))}
Nous allons tester nos régressions sur des processus dont nous connaissons les paramètres. Cette approche nous permettra d’évaluer la qualité et la fonctionnalité de nos méthodes.
Pour effectuer une régression de Poisson homogène, la fonction bru se base sur une grille arbitraire couvrant l’ensemble du domaine. Elle prend comme argument un SpatialPointsDataFrame contenant (pour coordonnées) les centres des cellules de cette grille et (pour valeur) le nombre de points du semis contenus dans chaque cellule correspondante. Il nous faut donc au préalable construire cette grille et son DataFrame.
# Simulation
poish = preparation(rpoispp(lambda = 10, win=bnd_owin), mesh = FALSE)
# Grille pour la régression de poisson
nb_carres_par_cote = 10 # résolution de la grille (arbitraire)
grid = raster(xmn=0, ymn=0, xmx=5, ymx=5, res=5/nb_carres_par_cote) # création de la grille sur le domaine
grid[] = 0 # grid[] : nombre d'observations dans chaque cellule
tab = table(cellFromXY(grid, poish$coord_pts)) # compte pour chaque cellule de la grille, le nombre d'observations
grid[as.numeric(names(tab))] = tab # qu'elle contient,
count_df = data.frame(coordinates(grid), count=grid[]) # ce résultat est ensuite stocké dans la variable count
coordinates(count_df) <- c("x", "y") # conversion en SpatialPointsDataFrame
# Régression
fit_poish = bru(count ~ Intercept, count_df, family="poisson")
lambda_poish = predict(fit_poish, pix, ~ exp(Intercept))
lambda_poish_df = as.data.frame(lambda_poish)
# Pour chaque point de pix (en ligne), on a la loi du lambda associé.
lambda_poish_df$intensity = lambda_poish$mean*(nb_carres_par_cote**2/25)
# Renormalisation pour avoir l'intensité, au lieu du nombre de points par cellule de la grille.
# Plot
ggplot() + gg(poish$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
plot(grid); points(poish$coord_pts, pch=20)
# Valeurs
print(paste("Intensité théorique 10"))
## [1] "Intensité théorique 10"
print(paste("Intensité observée", poish$nb_obs/25))
## [1] "Intensité observée 9.76"
print(paste("Intensité poisson homogène", exp(fit_poish$summary.fixed$mean)*(nb_carres_par_cote**2/25) ))
## [1] "Intensité poisson homogène 9.75994559179058"
summary(fit_poish)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): 3.899e+02
## Deviance Information Criterion (DIC): 3.896e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept 0.8919925 0.06407409 0.7638826 0.8927921 1.015553 0.8943809
## signif
## Intercept TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## None.
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## NULL
Nous souhaitons générer un processus non-homogène sur notre domaine. Poour cela nous choisissons arbitrairement comme intensité théorique sur le champ, la fonction qui somme les coordonées.
# Intensité théorique
intensity_pois = function(x,y){x+y}
grille = expand.grid(x=seq(0,5,length=200), y=seq(0,5,length=200))
grille$intensity = grille$x + grille$y
# Simulation
pois = preparation(rpoispp(lambda = intensity_pois, win=bnd_owin))
# Intensité observée
densi = density(pois$coord_pts_ppp, sigma=bw.diggle(pois$coord_pts_ppp))
# Echelle de couleurs commune
intensity_max <- ceiling(max(max(densi), max(grille$intensity)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)
# Plot
ggplot() + gg(pois$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
mat_pois = matrix(grille$intensity, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_pois, xcol=seq(ncol(mat_pois)), yrow=seq(nrow(mat_pois))), main="Intensité theorique", col=col_fixe)
plot(densi, main="Intensité observée", col=col_fixe)
Nous simulons ci-dessous un processus de Cox homogène muni d’une fonction de covariance de Matern, dont les paramètres sont \(moy. = 10\), \(var. = \exp(0.2)\), \(\alpha=1/2\) et \(\nu=1\).
# Simulation
lgcph = rLGCP(model='matern', mu = log(10), var=0.2, scale=1/2, nu=1, win=bnd_owin)
coxh = preparation(data.frame(x = lgcph$y, y = lgcph$x))
# lambda_theo est l'intensité exacte de notre processus et c'est une réalisation du champ aléatoire théorique de moyenne 10
lambda_theoh <- attr(lgcph, 'Lambda')
lambda_theoh_df = data.frame(expand.grid(x = lambda_theoh$yrow, y = lambda_theoh$xcol), val=as.vector(lambda_theoh$v))
# Régression
fit_coxh = lgcp(coordinates ~ Intercept, coxh$coord_pts, samplers = bnd)
# Plot
ggplot() + gg(coxh$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
ggplot(lambda_theoh_df, aes(x=x, y=y)) + geom_raster(aes(fill = val, interpolate=TRUE)) + coord_fixed() + colsc(max(lambda_theoh_df$val)) + gg(coxh$coord_pts) + ggtitle(label = "Lambda théorique")
## Warning: Ignoring unknown aesthetics: interpolate
# Valeurs
print(paste("Intensité théorique 10"))
## [1] "Intensité théorique 10"
print(paste("Intensité observée", coxh$nb_obs/25))
## [1] "Intensité observée 9.64"
print(paste("Intensité cox homogène", exp(fit_coxh$summary.fixed$mean)))
## [1] "Intensité cox homogène 9.70222535866406"
summary(fit_coxh)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -6.103e+02
## Deviance Information Criterion (DIC): -6.113e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept 2.272355 0.0644722 2.143435 2.273165 2.39667 2.274774
## signif
## Intercept TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## None.
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## NULL
Cette fois, nous simulons un processus de Cox non-homogène dont la moyenne n’est alors plus constante (nous choisissons la fonction damier de densité 5 et 20). Le processus est toujours muni d’une fonction de covariance de Matern dont les paramètres sont \(var. = \exp(0.2)\), \(\alpha=1/2\) et \(\nu=1\).
Afin d’effectuer la régression de Cox non-homogène, nous devons définir une variable de type champ aléatoire qui représentera les interactions intra-spécifique et nous permettra d’expliquer l’intensité. Le package INLA apporte une solution inattendue pour cette variable ; ce champ peut en effet être nommé comme on le souhaite et est reconnu seulement grâce à ces arguments (map = …., model = ….). Cette fonction est appelée “f” est dans la documentation R. Dans les exemples ci-dessous elle sera appelée “field”. Nous définissons “field” comme un champ aléatoire de covariance de Matérn, nous définissons donc un sous-modèle avec les paramètres du champ comme inconnus grâce à la fonction inla.spde2.pcmatern.
# Intensité exacte
fct_damier = function(x,y){n = length(x); return(log(5+15*(mod(floor(x)+floor(y), rep(2,n)))))}
grille = expand.grid(x=seq(0,5,length=200), y=seq(0,5,length=200))
grille$intensity = exp(fct_damier(grille$x, grille$y))
# Simulation
lgcp_sim = rLGCP(model='matern', mu=fct_damier, var=0.2, scale=1/2, nu=1, win=bnd_owin)
cox = preparation(data.frame(x = lgcp_sim$y, y = lgcp_sim$x))
# lambda_theo est l'intensité exacte de notre processus et c'est une réalisation du champ aléatoire théorique de moyenne 10
lambda_theo <- attr(lgcp_sim, 'Lambda')
lambda_theo_df = data.frame(expand.grid(x = lambda_theo$yrow, y = lambda_theo$xcol), val=as.vector(lambda_theo$v))
# Régression
model_matern = inla.spde2.pcmatern(cox$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, cox$coord_pts, samplers = bnd)
lambda_cox = predict(fit_cox, pix, ~ exp(field + Intercept))
lambda_cox_df = as.data.frame(lambda_cox)
# Plot
ggplot() + gg(cox$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
ggplot(grille, aes(x=x, y=y)) + geom_raster(aes(fill = intensity, interpolate=TRUE)) + coord_fixed() + colsc(max = max(lambda_theo_df$val)) + ggtitle(label = "Intensité théorique")
## Warning: Ignoring unknown aesthetics: interpolate
ggplot(lambda_theo_df, aes(x=x, y=y)) + geom_raster(aes(fill = val, interpolate=TRUE)) + coord_fixed() + colsc(max = max(lambda_theo_df$val)) + gg(cox$coord_pts) + ggtitle(label = "Intensité de la réalisation du champ")
## Warning: Ignoring unknown aesthetics: interpolate
plot(density(cox$coord_pts_ppp, sigma = bw.diggle(cox$coord_pts_ppp)), main="Intensité observée")
mat_cox = matrix(lambda_cox_df$mean, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_cox, xcol=seq(ncol(mat_cox)), yrow=seq(nrow(mat_cox))), main="Intensité moyenne")
mat_cox = matrix(lambda_cox_df$sd, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_cox, xcol=seq(ncol(mat_cox)), yrow=seq(nrow(mat_cox))), main="Ecart-type")
summary(fit_cox)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -8.989e+02
## Deviance Information Criterion (DIC): -8.991e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept 2.489473 0.5773392 1.381838 2.489828 3.588412 2.492595
## signif
## Intercept TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-0.00851155, 0.00603428], sd = [0.526497, 0.563169], quantiles = [-1.06471 : 1.06777]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 42.32952059 58.69558380 4.020955550 24.89737502 188.400863
## Stdev for field 0.03129879 0.03642442 0.002680121 0.02035396 0.126565
## mode
## Range for field 10.179032859
## Stdev for field 0.007305385
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq
## 1 range 41.848345415 2.871187e+03 53.583455583 3.893303e+00
## 2 log.range 3.244986768 9.562635e-01 0.977887256 1.399332e+00
## 3 variance 0.002147927 4.871517e-05 0.006979626 -6.087022e-05
## 4 log.variance -7.848836246 3.848233e+00 1.961691335 -1.183507e+01
## median uq
## 1 24.6982490732 185.14378766
## 2 3.2115729110 5.22760208
## 3 0.0003483903 0.01525113
## 4 -7.7952687180 -4.15498732
La fonction qui simule un processus de Matérn est rMatClust(kappa, scale, mu, win, …. ) de spatstat avec :
kappa = intensité du premier processus de poisson (celui qui génère les centres des cercles) scale = rayon des cercles mu = nombre de points moyen par cercle win = fenêtre globale
# Simulation
matern = rMatClust(kappa=0.3, scale=0.4, mu=20, win = bnd_owin)
mat = preparation(data.frame(x = matern$x, y = matern$y))
# Régression
model_matern = inla.spde2.pcmatern(mat$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_mat = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, mat$coord_pts, samplers = bnd)
lambda_mat = predict(fit_mat, pix, ~ exp(field + Intercept))
lambda_mat_df = as.data.frame(lambda_mat)
# Plot
ggplot() + gg(mat$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
plot(density(mat$coord_pts_ppp, sigma = bw.diggle(mat$coord_pts_ppp)), main="Intensité observée")
mat_mat = matrix(lambda_mat_df$mean, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_mat, xcol=seq(ncol(mat_mat)), yrow=seq(nrow(mat_mat))), main="Intensité moyenne")
mat_sd = matrix(lambda_mat_df$sd, ncol=200, nrow=200, byrow = TRUE)
plot(im(mat_sd, xcol=seq(ncol(mat_sd)), yrow=seq(nrow(mat_sd))), main="Ecart-type")
summary(fit_mat)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -9.579e+02
## Deviance Information Criterion (DIC): -9.598e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept 1.439986 0.4833149 0.4907232 1.440095 2.387725 1.440355
## signif
## Intercept TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-2.07012, 3.11968], sd = [0.553607, 1.6963], quantiles = [-5.21758 : 4.46431]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 40.4884204 60.4016291 3.106282508 22.66336008 188.5357188
## Stdev for field 0.0418197 0.0620962 0.002669162 0.02331727 0.1959196
## mode
## Range for field 8.041174555
## Stdev for field 0.007020902
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq
## 1 range 39.969590947 2.982100e+03 54.60860853 2.955772819
## 2 log.range 3.140773192 1.087530e+00 1.04284728 1.141436642
## 3 variance 0.004894685 4.806958e-04 0.02192478 -0.007111195
## 4 log.variance -7.528669521 4.758744e+00 2.18145466 -11.837313864
## median uq
## 1 22.457889591 185.0392461
## 2 3.117215524 5.2268952
## 3 0.000016481 0.0274249
## 4 -7.523275155 -3.2879265
Les données d’altitude dont nous disposons sont les altitudes au pied de chaque arbre de la parcelle.
# Altitude pour tous les arbres
altitude = data.frame(x=p16$X,y=p16$Y,z=p16$altitude)
coordinates(altitude)=c("x","y")
ggplot(p16, aes(x=X, y=Y)) + geom_point(aes(color = altitude)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle("Altitude exacte pour tous les arbres")
# Altitude sur un échantillon d'arbres dans la parcelle (nécessaire à la rapidité du krigeage)
diviseur = 10
sample_p16 = p16[sample(1:nrow(p16), floor(nrow(p16)/diviseur)),]
ggplot(sample_p16, aes(x=X, y=Y)) + geom_point(aes(color = altitude)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle("Altitude exacte sur un échantillon des arbres")
# Visualisation sur la carte
library(rgdal)
## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/lena.klay/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/lena.klay/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1
library(leaflet)
library(mapview)
coord_utm = data.frame(long=p16$Xutm, lat=p16$Yutm)
utms = SpatialPoints(coord_utm, proj4string=CRS('+proj=utm +zone=22'))
longlats <- spTransform(utms, CRS('+init=epsg:4326'))
longlats$alt = p16$altitude
mapview(longlats, zcol = c("alt"), legend = TRUE)
(Attention ici : passage du carré de 500x500 au carré de 5x5)
Nous souhaitons rajouter l’altitude comme variable du modèle ; non pas comme un champ aléatoire mais comme un champ fixé. La syntaxe est donc légèrement différente, la fonction “f” (qu’on nommera ici “alti” comme altitude) prendra en argument model = “linear”.
Toujours dans cette même fonction, l’argument map est cette fois défini comme une seconde fonction qui retourne l’altitude à chaque coordonnée des points considérés. Cette seconde fonction nécessite un SpatialPixelsDataFrame de l’altitude, c’est à dire des valeurs d’altitude placées sur une grille régulière, ce qui n’est actuellement pas le cas puisque nos données correspondent à l’emplacement de chaque arbre. Nous devons donc effectuer au préalable un krigeage, détaillé ci-dessous.
Attention : si vous obtenez une erreur de type “Error in match.names(clabs, names(xi)) : names do not match previous names”, essayez de changer le nom de la fonction “f”… tous les noms ne sont pas acceptés et les messages d’erreur ne sont pas très clairs.
# Préparation des données
altitude = SpatialPointsDataFrame(coords = data.frame(x=p16$X/100,y=p16$Y/100),
data = data.frame(z=p16$altitude), proj4string = CRS(as.character(NA)))
sample_altitude = SpatialPointsDataFrame(coords = data.frame(x=sample_p16$X/100,y=sample_p16$Y/100),
data = data.frame(z=sample_p16$altitude), proj4string = CRS(as.character(NA)))
sample_altitude <- sample_altitude[-zerodist(sample_altitude)[,1],] # enlève les points de même coordonnées
# Krigeage
vgmEmpirique <- gstat::variogram(z ~ 1, data = sample_altitude)
vgmX <- fit.variogram(vgmEmpirique, vgm("Gau")) # Ajustement d'un modèle gaussien
geoX <- gstat(formula = z ~ 1, locations = sample_altitude, model = vgmX) # Objet geostat qui décrit toutes les # caractéristiques de la modélisation.
grid <- expand.grid((0:100)/20, (0:100)/20) # Krigeage sur une grille de résolution 5cm (carré de 5m x 5m).
names(grid) <- c("x", "y")
gridded(grid) <- ~x + y
altitude_grid <- predict(geoX, newdata = grid) # Calcul de la valeur de l'altitude sur les points de la grille
## [using ordinary kriging]
altitude_grid=altitude_grid[,-2]
names(altitude_grid) = c("z")
# Plot
image(altitude_grid, col = topo.colors(20, alpha = 1), asp = 1)
contour(altitude_grid, add = TRUE)
title(main = "Krigeage de l'altitude", font.main = 4)
Seconde fonction qui retourne l’altitude des points étudiés
fct_alt <- function(x,y) { # /!\ altitude_grid DOIT ABSOLUMENT être un SpatialPixelsDataFrame
spp <- SpatialPoints(data.frame(x=x,y=y)) # SpatialPoint object
proj4string(spp) <- CRS(proj4string(altitude_grid)) # auquel on attache le système de coordonnées de l'objet altitude_grid
v <- over(spp, altitude_grid) # extrait l'altitude correspondante aux coordonnées des points
v[is.na(v)] <- 0 # enlève les NA
return(v$z)
}
Nous testons ensuite notre modèle comprenant l’altitude sur quelques espèces d’arbres.
# Choix de l'espèce et du genre
genre = "Vouacapoua"
espece = "americana"
# Préparation des données
voam = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))
# Altitude seule
fit_alt_voam = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = voam$coord_pts, samplers = bnd)
lambda_alt_voam <- predict(fit_alt_voam, pix, ~ exp(alti + Intercept))
# Champ SPDE seul
model_matern = inla.spde2.pcmatern(voam$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_voam = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = voam$coord_pts, samplers = bnd)
lambda_spde_voam <- predict(fit_spde_voam, pix, ~ exp(field + Intercept))
# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(voam$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_voam = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = voam$coord_pts, samplers = bnd)
lambda_alt_spde_voam <- predict(fit_alt_spde_voam, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_voam_df = as.data.frame(lambda_alt_spde_voam)
# Plot
# Semis de points
ggplot() + gg(voam$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))
# Intensité moyenne prédite des 3 modèles
max_voam = max(max(lambda_alt_voam$mean),max(lambda_spde_voam$mean),max(lambda_alt_spde_voam$mean)) # Max. échelle de couleurs
ggplot() + gg(lambda_alt_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_voam)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_spde_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_voam)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_alt_spde_voam) + gg(bnd) + gg(voam$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_voam)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_voam, "alti")
plot(fit_alt_spde_voam, "alti")
# Critères DIC et WAIC
scores <- data.frame(
dic=c(alt = fit_alt_voam$waic$waic, spde = fit_spde_voam$waic$waic, alt_spde = fit_alt_spde_voam$waic$waic),
waic=c(alt = fit_alt_voam$dic$dic, spde = fit_spde_voam$dic$dic, alt_spde = fit_alt_spde_voam$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores
## dic waic
## altitude -512.9292 -515.3266
## field -566.1474 -566.5644
## altitude + field -570.1489 -570.5601
summary(fit_alt_voam)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -5.129e+02
## Deviance Information Criterion (DIC): -5.153e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti 0.1982447 0.01987555 0.1597474 0.1980681 0.237728367
## Intercept -0.6415450 0.33673873 -1.3205301 -0.6352843 0.002163823
## mode signif
## alti 0.1977187 TRUE
## Intercept -0.6228578 FALSE
##
## --- Random effects -------------------------------------------------------------------------------
##
## None.
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## NULL
summary(fit_spde_voam)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -5.661e+02
## Deviance Information Criterion (DIC): -5.666e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept 0.8946979 0.6157455 -0.3557254 0.9033231 2.099715 0.9210275
## signif
## Intercept FALSE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-1.82194, 3.13234], sd = [0.641126, 1.05395], quantiles = [-3.96713 : 4.55821]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 4.4804277 0.91161583 2.9475646 4.3932478 6.5215763
## Stdev for field 0.5897562 0.08069693 0.4479424 0.5840389 0.7643381
## mode
## Range for field 4.2259591
## Stdev for field 0.5726136
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq median
## 1 range 4.4792521 0.818602944 0.9047668 2.9552339 4.3917884
## 2 log.range 1.4793818 0.040677686 0.2016871 1.0819332 1.4794914
## 3 variance 0.3541593 0.009519242 0.0975666 0.2010395 0.3409422
## 4 log.variance -1.0746004 0.074010033 0.2720479 -1.6062641 -1.0763325
## uq
## 1 6.4953877
## 2 1.8723103
## 3 0.5817151
## 4 -0.5400997
summary(fit_alt_spde_voam)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -5.701e+02
## Deviance Information Criterion (DIC): -5.706e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti 0.1152663 0.02337132 0.06936644 0.1152541 0.1611970
## Intercept -0.5149656 0.67676818 -1.86962484 -0.5108825 0.8168936
## mode signif
## alti 0.1152328 TRUE
## Intercept -0.5020177 FALSE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-1.7757, 2.37553], sd = [0.620763, 0.894035], quantiles = [-3.61457 : 3.76245]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 5.2571960 1.168258 3.330659 5.1333402 7.9050183
## Stdev for field 0.4874612 0.075589 0.355482 0.4817944 0.6522584
## mode
## Range for field 4.8961122
## Stdev for field 0.4708139
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq median
## 1 range 5.2555776 1.343533145 1.15910877 3.3382363 5.1314740
## 2 log.range 1.6354701 0.048247758 0.21965372 1.2036900 1.6351308
## 3 variance 0.2431839 0.005784601 0.07605656 0.1268349 0.2319968
## 4 log.variance -1.4608678 0.095105293 0.30839146 -2.0670681 -1.4613436
## uq
## 1 7.8707730
## 2 2.0644874
## 3 0.4230275
## 4 -0.8584636
Etude approfondie du modèle field + altitude dans le cas de Vouacapoua americana :
# Intensité observée avec density
densi = density(voam$coord_pts_ppp, sigma = bw.diggle(voam$coord_pts_ppp))
densi$v[densi$v<0]=0
# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_voam_df$mean, ncol=200, nrow=200, byrow = TRUE)
# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_voam_df$sd, ncol=200, nrow=200, byrow = TRUE)
# Echelle de couleurs commune pour intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_voam_df$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)
# Plot
plot(densi, main="Intensité observée", col=col_fixe)
plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)
plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type")
post.range = spde.posterior(fit_alt_spde_voam, name="field", what="range"); plot(post.range)
post.matcorr = spde.posterior(fit_alt_spde_voam, name="field",what="matern.correlation"); plot(post.matcorr)
# Choix de l'espèce et du genre
genre = "Eperua"
espece = "falcata"
# Préparation des données
epfa = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))
# Altitude seule
fit_alt_epfa = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = epfa$coord_pts, samplers = bnd)
lambda_alt_epfa <- predict(fit_alt_epfa, pix, ~ exp(alti + Intercept))
# Champ SPDE seul
model_matern = inla.spde2.pcmatern(epfa$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_epfa = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = epfa$coord_pts, samplers = bnd)
lambda_spde_epfa <- predict(fit_spde_epfa, pix, ~ exp(field + Intercept))
# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(epfa$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_epfa = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = epfa$coord_pts, samplers = bnd)
lambda_alt_spde_epfa <- predict(fit_alt_spde_epfa, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_df_epfa = as.data.frame(lambda_alt_spde_epfa)
# Plot
# Semis de points
ggplot() + gg(epfa$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))
# Intensité moyenne prédite des 3 modèles
max_epfa = max(max(lambda_alt_epfa$mean),max(lambda_spde_epfa$mean),max(lambda_alt_spde_epfa$mean)) # max. échelle de couleurs
ggplot() + gg(lambda_alt_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_epfa)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_spde_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_epfa)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_alt_spde_epfa) + gg(bnd) + gg(epfa$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_epfa)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_epfa, "alti")
plot(fit_alt_spde_epfa, "alti")
# Critères DIC et WAIC
scores <- data.frame(
dic=c(alt = fit_alt_epfa$waic$waic, spde = fit_spde_epfa$waic$waic, alt_spde = fit_alt_spde_epfa$waic$waic),
waic=c(alt = fit_alt_epfa$dic$dic, spde = fit_spde_epfa$dic$dic, alt_spde = fit_alt_spde_epfa$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores
## dic waic
## altitude -143.1317 -143.7463
## field -609.1629 -607.6059
## altitude + field -608.5919 -606.8954
summary(fit_alt_epfa)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -1.431e+02
## Deviance Information Criterion (DIC): -1.437e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti -0.06904678 0.02043147 -0.1091926 -0.0690363 -0.02899805
## Intercept 2.41347156 0.23360976 1.9413661 2.4182071 2.85894137
## mode signif
## alti -0.06901356 TRUE
## Intercept 2.42760840 TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## None.
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## NULL
summary(fit_spde_epfa)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -6.092e+02
## Deviance Information Criterion (DIC): -6.076e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Intercept -0.3502781 0.6574403 -1.697013 -0.3326241 0.899974 -0.2975257
## signif
## Intercept FALSE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-1.77624, 5.74104], sd = [0.717743, 2.02388], quantiles = [-5.4864 : 8.35587]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field 2.315757 0.3129724 1.7667776 2.293108 2.992918 2.247193
## Stdev for field 1.142503 0.1078816 0.9458953 1.137194 1.369543 1.126442
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq median
## 1 range 2.3155014 0.09664505 0.3108779 1.7686429 2.2926387
## 2 log.range 0.8307180 0.01802728 0.1342657 0.5690927 0.8295366
## 3 variance 1.3166633 0.06183539 0.2486672 0.8963457 1.2928225
## 4 log.variance 0.2575794 0.03545913 0.1883059 -0.1109539 0.2566033
## uq
## 1 2.9874274
## 2 1.0952718
## 3 1.8698500
## 4 0.6270328
summary(fit_alt_spde_epfa)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -6.086e+02
## Deviance Information Criterion (DIC): -6.069e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti -0.03269948 0.05442809 -0.1408992 -0.03231621 0.07325973
## Intercept 0.03650594 0.91217420 -1.8273892 0.06175723 1.75923485
## mode signif
## alti -0.03156701 FALSE
## Intercept 0.11160607 FALSE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-1.73656, 5.70725], sd = [0.721576, 2.03381], quantiles = [-5.43981 : 8.25936]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for field 2.274916 0.3138912 1.7250694 2.252120 2.954831 2.206272
## Stdev for field 1.143510 0.1087392 0.9440265 1.138576 1.371684 1.129071
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq median
## 1 range 2.274655 0.09721038 0.3117858 1.7266341 2.2516470
## 2 log.range 0.812541 0.01879556 0.1370969 0.5450284 0.8114918
## 3 variance 1.319140 0.06284580 0.2506907 0.8934717 1.2959568
## 4 log.variance 0.259194 0.03602727 0.1898085 -0.1141936 0.2590217
## uq
## 1 2.9489123
## 2 1.0823099
## 3 1.8748176
## 4 0.6296763
Etude approfondie du modèle field + altitude dans le cas de Eperua falcata :
# Intensité observée avec density
densi = density(epfa$coord_pts_ppp, sigma = bw.diggle(epfa$coord_pts_ppp))
densi$v[densi$v<0]=0
# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_df_epfa$mean, ncol=200, nrow=200, byrow = TRUE)
# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_df_epfa$sd, ncol=200, nrow=200, byrow = TRUE)
# Echelle de couleurs commune pour l'intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_df_epfa$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)
# Plot
plot(densi, main="Intensité observée", col=col_fixe)
plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)
plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type")
post.range = spde.posterior(fit_alt_spde_epfa, name="field", what="range"); plot(post.range)
post.matcorr = spde.posterior(fit_alt_spde_epfa, name="field",what="matern.correlation"); plot(post.matcorr)
# Choix de l'espèce et du genre
genre = "Oenocarpus"
espece = "bataua"
# Préparation des données
oeba = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece]/100, y = p16$Y[p16$Genus==genre&p16$Species==espece]/100))
# Altitude seule
fit_alt_oeba = lgcp(coordinates ~ alti(map = fct_alt(x,y), model = "linear") + Intercept, data = oeba$coord_pts, samplers = bnd)
lambda_alt_oeba <- predict(fit_alt_oeba, pix, ~ exp(alti + Intercept))
# Champ SPDE seul
model_matern = inla.spde2.pcmatern(oeba$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_spde_oeba = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = oeba$coord_pts, samplers = bnd)
lambda_spde_oeba <- predict(fit_spde_oeba, pix, ~ exp(field + Intercept))
# Altitude + champ SPDE
model_matern = inla.spde2.pcmatern(oeba$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_alt_spde_oeba = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + alti(map = fct_alt(x,y), model = "linear") + Intercept, data = oeba$coord_pts, samplers = bnd)
lambda_alt_spde_oeba <- predict(fit_alt_spde_oeba, pix, ~ exp(field + alti + Intercept))
lambda_alt_spde_df_oeba = as.data.frame(lambda_alt_spde_oeba)
# Plot
# Semis de points
ggplot() + gg(oeba$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
# Points avec altitude
abc = data.frame(a=p16$X[p16$Genus==genre&p16$Species==espece]/100,b=p16$Y[p16$Genus==genre&p16$Species==espece]/100,c=p16$altitude[p16$Genus==genre&p16$Species==espece])
ggplot(abc, aes(x=a, y=b)) + geom_point(aes(color = c)) + scale_color_viridis(option="D") + coord_fixed() + ggtitle(paste(genre, espece))
# Intensité moyenne prédite des 3 modèles
max_oeba = max(max(lambda_alt_oeba$mean),max(lambda_spde_oeba$mean),max(lambda_alt_spde_oeba$mean)) # max. échelle de couleurs
ggplot() + gg(lambda_alt_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude") + colsc(max_oeba)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_spde_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("SPDE") + colsc(max_oeba)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
ggplot() + gg(lambda_alt_spde_oeba) + gg(bnd) + gg(oeba$coord_pts, shape="+") + coord_equal() + ggtitle("Altitude + SPDE") + colsc(max_oeba)
## Regions defined for each Polygons
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
# Loi à posteriori des coefficients linéaires de l'altitude
plot(fit_alt_oeba, "alti")
plot(fit_alt_spde_oeba, "alti")
# Critères DIC et WAIC
scores <- data.frame(
dic=c(alt = fit_alt_oeba$waic$waic, spde = fit_spde_oeba$waic$waic, alt_spde = fit_alt_spde_oeba$waic$waic),
waic=c(alt = fit_alt_oeba$dic$dic, spde = fit_spde_oeba$dic$dic, alt_spde = fit_alt_spde_oeba$dic$dic))
rownames(scores) <- c("altitude", "field", "altitude + field")
scores
## dic waic
## altitude -901.5785 -903.5977
## field -933.0700 -933.6494
## altitude + field -933.2780 -933.7714
summary(fit_alt_oeba)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -9.016e+02
## Deviance Information Criterion (DIC): -9.036e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti -0.01179511 0.01206122 -0.03544741 -0.01180521 0.01189162
## Intercept 2.64143243 0.15674084 2.32830073 2.64329762 2.94393736
## mode signif
## alti -0.01182442 FALSE
## Intercept 2.64700462 TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## None.
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## NULL
summary(fit_spde_oeba)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -9.331e+02
## Deviance Information Criterion (DIC): -9.336e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant mode signif
## Intercept 2.459096 0.3387078 1.759711 2.461535 3.150623 2.46739 TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-0.528934, 0.491626], sd = [0.3415, 0.411035], quantiles = [-1.40309 : 1.30928]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 7.9403136 2.6526225 4.08837086 7.4846036 14.3944840
## Stdev for field 0.1767393 0.0477657 0.09970707 0.1712971 0.2861291
## mode
## Range for field 6.6663295
## Stdev for field 0.1609388
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq
## 1 range 7.93466725 6.8838889509 2.62371663 4.106581538
## 2 log.range 2.02025022 0.1022367481 0.31974482 1.410394843
## 3 variance 0.03345191 0.0003457436 0.01859418 0.009957444
## 4 log.variance -3.53770321 0.2877557107 0.53642866 -4.611511038
## median uq
## 1 7.47955216 14.30889014
## 2 2.01173690 2.66281840
## 3 0.02930361 0.08105774
## 4 -3.53050081 -2.50984618
summary(fit_alt_spde_oeba)
##
## --- Likelihoods ----------------------------------------------------------------------------------
##
##
## --- Criteria -------------------------------------------------------------------------------------
##
## Watanabe-Akaike information criterion (WAIC): -9.333e+02
## Deviance Information Criterion (DIC): -9.338e+02
##
## --- Fixed effects --------------------------------------------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## alti -0.02434751 0.01451099 -0.05299715 -0.02430024 0.004005865
## Intercept 2.76132361 0.37867946 1.99356014 2.76410644 3.514588434
## mode signif
## alti -0.02420687 FALSE
## Intercept 2.76985799 TRUE
##
## --- Random effects -------------------------------------------------------------------------------
##
## field ranges: mean = [-0.558687, 0.423791], sd = [0.337832, 0.415136], quantiles = [-1.44352 : 1.22923]
##
## --- All hyper parameters (internal representation) -----------------------------------------------
##
## mean sd 0.025quant 0.5quant 0.975quant
## Range for field 7.9407660 2.67227909 4.06742403 7.4796476 14.4480425
## Stdev for field 0.1745113 0.04814275 0.09698259 0.1689933 0.2847932
## mode
## Range for field 6.6529763
## Stdev for field 0.1584815
##
##
## --- Field 'field' transformed hyper parameters ---
## param mean var sd lq median
## 1 range 7.93503961 6.9853394734 2.64297928 4.08581952 7.47456187
## 2 log.range 2.01958195 0.1036717215 0.32198093 1.40531969 2.01106915
## 3 variance 0.03270457 0.0003445674 0.01856253 0.00943151 0.02851929
## 4 log.variance -3.56611774 0.3003139629 0.54800909 -4.66568056 -3.55762651
## uq
## 1 14.36165942
## 2 2.66651263
## 3 0.08033236
## 4 -2.51880418
Etude approfondie du modèle field + altitude dans le cas de Oenocarpus bataua :
# Intensité observée avec density
densi = density(oeba$coord_pts_ppp, sigma = bw.diggle(oeba$coord_pts_ppp))
densi$v[densi$v<0]=0
# Intensité moyenne prédite
mat_int = matrix(lambda_alt_spde_df_oeba$mean, ncol=200, nrow=200, byrow = TRUE)
# Ecart-type de la prédiction
mat_eca = matrix(lambda_alt_spde_df_oeba$sd, ncol=200, nrow=200, byrow = TRUE)
# Echelle de couleurs commune pour l'intensité obs. et moy. prédite
intensity_max <- ceiling(max(max(densi), max(lambda_alt_spde_df_oeba$mean)))
col_fixe <- colourmap(hcl.colors(intensity_max), breaks=0:intensity_max)
# Plot
plot(densi, main="Intensité observée", col=col_fixe)
plot(im(mat_int, xcol=seq(ncol(mat_int)), yrow=seq(nrow(mat_int))), main="Intensité moy. prédite", col=col_fixe)
plot(im(mat_eca, xcol=seq(ncol(mat_eca)), yrow=seq(nrow(mat_eca))), main="Ecart-type")
post.range = spde.posterior(fit_alt_spde_oeba, name="field", what="range"); plot(post.range)
post.matcorr = spde.posterior(fit_alt_spde_oeba, name="field",what="matern.correlation"); plot(post.matcorr)
Comme expliqué dans le rapport, la méthode INLA-SPDE ne permet pas d’étudier les interactions entre les espèces sur la parcelle. Des outils non-paramètriques plus classiques nous donne un aperçu de ce phénomène.
# Préparation des deux espèces
genre = "Qualea"
espece = "rosea"
quro_grd = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece], y = p16$Y[p16$Genus==genre&p16$Species==espece]), size=500)
genre = "Vouacapoua"
espece = "americana"
voam_grd = preparation(data.frame(x = p16$X[p16$Genus==genre&p16$Species==espece], y = p16$Y[p16$Genus==genre&p16$Species==espece]), size=500)
# Graphique des deux espèces
doubleplot = data.frame( rbind(quro_grd$df,voam_grd$df) , Espece = c(rep("Qualea rosea",quro_grd$nb_obs),rep("Vouacapoua americana",voam_grd$nb_obs)) )
ggplot(doubleplot, aes(x=x, y=y)) + geom_point(aes(color = Espece)) + coord_fixed()
# Calcul de la fonction M
X = c(voam_grd$df[,1],quro_grd$df[,1])
Y = c(voam_grd$df[,2],quro_grd$df[,2])
PointType <- c(rep("V. Americana", voam_grd$nb_obs), rep("Q. Rosea", quro_grd$nb_obs))
PointWeight <- rep(1,voam_grd$nb_obs+quro_grd$nb_obs)
pattern <- wmppp(data.frame(X, Y, PointType, PointWeight), owin(xrange=c(0,500), yrange=c(0,500)))
result_M = Mhat(X=pattern, ReferenceType="V. Americana", NeighborType="Q. Rosea")
plot(result_M$r, result_M$M, type = "l", ylim=c(0.15,1), xlab="r", ylab="M(r)")
lines(result_M$r, result_M$theo, type = "l", col = "red", lty=2)
legend(x=250, y=0.32, legend=c("M indépendant", "M observé"), col=c("red", "black"), lty=c(2,1))
Comme détaillé dans le rapport, la fonction lgcp présente une incohérence dans son implémentation. Afin de mettre celle-ci en évidence, nous testons la fonction sur des processus bien pensés.
# Processus n°1 : Poisson homogène (lambda=20) sur la moitié du domaine
poish_moitie = rpoispp(lambda = 20, win = owin(xrange=c(0,2.5), yrange=c(0,5)))
processus_1 = preparation(data.frame(x=poish_moitie$x, y= poish_moitie$y))
# Processus n°2 : processus n°1 auquel on a rajouté deux points aux deux coins opposés du carré (de coord. (5,0) et (5,5)).
poish_moitie_deux_coins = data.frame(x = c(processus_1$df$x, 5, 5), y= c(processus_1$df$y, 0, 5))
processus_2 = preparation(data.frame(x=poish_moitie_deux_coins$x, y= poish_moitie_deux_coins$y))
# Processus n°3 : processus n°2 mais dont les points ne sont cette fois plus sur la frontière mais légèrement rentrés dans le domaine (de coord. (4.8,0.2) et (4.8,4.8)).
poish_moitie_deux_coins_bis = data.frame(x = c(processus_1$df$x, 4.8, 4.8), y= c(processus_1$df$y, 0.2, 4.8))
processus_3 = preparation(data.frame(x=poish_moitie_deux_coins_bis$x, y= poish_moitie_deux_coins_bis$y))
# Régression de Cox homogène
fit_coxh1 <- lgcp(coordinates ~ Intercept, data = processus_1$coord_pts, samplers = bnd)
fit_coxh2 <- lgcp(coordinates ~ Intercept, data = processus_2$coord_pts, samplers = bnd)
fit_coxh3 <- lgcp(coordinates ~ Intercept, data = processus_2$coord_pts, samplers = bnd)
# Plot
ggplot() + gg(processus_1$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
ggplot() + gg(processus_2$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
ggplot() + gg(processus_3$coord_pts) + gg(bnd) + coord_fixed()
## Regions defined for each Polygons
# Intensités
print(paste("Intensité observée processus n°1 = ", (processus_1$nb_obs)/25))
## [1] "Intensité observée processus n°1 = 10.88"
print(paste("Intensité cox homogène processus n°1 = ", exp(fit_coxh1$summary.fixed$mean)))
## [1] "Intensité cox homogène processus n°1 = 19.8342735787074"
print(paste("Intensité observée processus n°2 = ", (processus_2$nb_obs)/25))
## [1] "Intensité observée processus n°2 = 10.96"
print(paste("Intensité cox homogène processus n°2 = ", exp(fit_coxh2$summary.fixed$mean)))
## [1] "Intensité cox homogène processus n°2 = 10.9596934073662"
print(paste("Intensité observée processus n°3 = ", (processus_3$nb_obs+2)/25))
## [1] "Intensité observée processus n°3 = 11.04"
print(paste("Intensité cox homogène processus n°3 = ", exp(fit_coxh3$summary.fixed$mean)))
## [1] "Intensité cox homogène processus n°3 = 10.9596934073662"
L’intensité cox homogène processus n°1 est généralement légèrement inférieure à 20 (intensité du poisson sur la moitié gauche du domaine) : la fonction ne semble pas prendre en considération les zones vides du domaine.
# Régression de Cox non-homogène processus n°1
model_matern = inla.spde2.pcmatern(processus_1$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox1 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_1$coord_pts, samplers = bnd)
lambda_cox1 <- predict(fit_cox1, pix, ~ exp(field + Intercept))
# Régression de Cox non-homogène processus n°2
model_matern = inla.spde2.pcmatern(processus_2$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox2 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_2$coord_pts, samplers = bnd)
lambda_cox2 <- predict(fit_cox2, pix, ~ exp(field + Intercept))
# Régression de Cox non-homogène processus n°3
model_matern = inla.spde2.pcmatern(processus_3$mesh, prior.sigma = c(0.1, 0.01), prior.range = c(5, 0.01))
fit_cox3 = lgcp(coordinates ~ field(map = coordinates, model = model_matern) + Intercept, data = processus_3$coord_pts, samplers = bnd)
lambda_cox3 <- predict(fit_cox3, pix, ~ exp(field + Intercept))
# Plot
plot(lambda_cox1)
plot(lambda_cox2)
plot(lambda_cox3)
Le problème d’implémentation n’impacte cependant pas la régression lorsque celle-ci est non-homogène.